package algorithm.sphericharmonics;

import algorithm.fft.Complex;

/*
 * The scaling of the polynomials is given here:

                   ________
                  / (l-m)!
     Q_{l,m} =   / -------- P_{l,m}
               \/   (l+m)!


    where P_{l,m} are the associated Legendre polynomials  and
    Q_{l,m}  are  the  scaled associated Legendre polynomials.
    The calculation of the polynomials is done with  the  fol-
    lowing three term recursion on l:

                               2  m/2
     Q_{m,m}(x) = c1(m)*( 1 - x  )

     Q_{m+1,m}(x) = c2(m)*x*Q_{m,m}(x)

     Q_{l,m}(x) =  c3(l,m)*x*Q_{l-1,m}(x)
                       - c4(l,m)*Q_{l-2,m}(x)

    where the functions of l and m are given below:

                        ____________________
                       /    m               |
                  m   /  --------'
     c1(m)  = (-1) * /  '  |  |          1
                    /      |  |    1 - -----
                   /       |  |         2*i
                 \/       i = 1

                 _____
                /     |
     c2(m) =   / 2*m+1
             \/

                 ( 2*l-1 )
     c3(l,m) = --------------
                   2   2  1/2
                ( l - m  )

                    _____________
                   /      2   2  |
                  /  (l-1) - m
     c4(l,m) =   / -------------
                /       2   2
              \/       l - m

 */
public final class ScaledAssociatedLagendrePolynomial {
	private int m; // should be always larger than or equal to 0; 
	public static final int lmax=VoxelizedSphere.lmax;
	private static final int[] lmIndex=VoxelizedSphere.lmIndex_offset; 
	private static final float[] sqrts_ll_mm=new float[((lmax+1)*lmax)/2];
	private static final float[] q_mm_factors=new float[lmax];
	static {
		for (int l=2; l<lmax;l++){
			for (int m=0;m<=l;m++){
				int index=lmIndex[l]+m;
				sqrts_ll_mm[index]=(float)Math.sqrt(l*l-m*m); 
			}
		}
		float q_mm_factor_m=1; 
		q_mm_factors[0]=1; 
		for (int m=1;m<lmax; m++){
			q_mm_factor_m=q_mm_factor_m*(float)Math.sqrt(1.0f-0.5f/m); 
			if (m%2!=0){
				q_mm_factors[m]=-q_mm_factor_m; 
			}else {
				q_mm_factors[m]=q_mm_factor_m;
			}
		}
	}
	private float x;
	protected float[] qValues; 
	public  ScaledAssociatedLagendrePolynomial(int rank,float val){
		m=rank;
		x=val;
		/*
		float v=1; 
		for (int i=1; i<=m; i++){
			v=v*(1.0f-1.0f/(2*i));  
		}
		v=(float)Math.sqrt(v); 
		if (m%2!=0){
			v=-v; 
		}
		if (Math.abs(v-q_mm_factors[m])>1.0e-5){
			throw new RuntimeException("wrong:"+"v/factor:"+v+"/"+q_mm_factors[m]+"/m:"+m); 
		}
		
		v=v*(float)Math.pow((1-x*x),((float)m)/2);
		*/
		
		qValues=new float[lmax+1-m]; // for the fixed m;
		qValues[0]=(float)q_mm_factors[m]*(float)Math.pow((1-x*x),((float)m)/2);
		qValues[1]=(float)Math.sqrt(2*m+1)*x*qValues[0]; 
	}
	public final float getPositiveValue(int l){
		return qValues[l-m]; 
	}
	public float[] getQValues(){
		return qValues; 
	}
	
	/*
	public static final float sqrt(float x){
		return (float)Math.sqrt(x); 
	}
	*/
	/*
	 * computer q(l=m,m),q(l=m+1,m), .... q(l=lmax,m); 
	 */
	public final void precompute(){
		for (int l=m+2; l<lmax; l++){
			float ql_1m=qValues[l-1-m];   //qlm(l-1,x);
			float ql_2m=qValues[l-2-m];   //qlm(l-2,x);
			int index=lmIndex[l]+m; 
			float sqrt_ll_mm=sqrts_ll_mm[index]; 
			//qValues[l-m]=(x*(2*l-1)*ql_1m-(float)Math.sqrt((l-1)*(l-1)-m*m)*ql_2m)/(float)Math.sqrt(l*l-m*m);
			qValues[l-m]=(x*(2*l-1)*ql_1m-(float)Math.sqrt((l-1)*(l-1)-m*m)*ql_2m)/sqrt_ll_mm;
		}
		/*
		for (int l=m; l<=lmax; l++){
			System.out.println(l+","+m+"):"+qValues[l-m]); 
		}
		*/
	}
	public static void main(String[] args){
		
		float x=(float)Math.PI/20; 
		int l=1; 
		ScaledAssociatedLagendrePolynomial y3_3=new ScaledAssociatedLagendrePolynomial(0,x);
		y3_3.precompute(); 
		System.out.println(y3_3.getPositiveValue(3));
		System.out.println(ScaledAssociatedLegenedrePolynomialFunct.pfunct(3,0,x));
	}
}
